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Section  1 


INTRODUCTION  AND  PROBLEM  DESCRIPTION 

Background 

In  1991  Dr.  John  Shoosmith,  a  researcher  at  the  National  Aeronautics  and 
Space  Administration  (NASA),  proposed  a  new  method  for  solving  systems  of 
first-order  differential  equations  in  parallel.  He  presented  his  method  and 
a  one-dimensional  application  of  the  algorithm  to  the  1991  International 
Conference  on  Industrial  and  Applied  Mathematics  held  in  Washington,  DC.  His 
original  research  proposal  indicated  plans  to  apply  the  method  to  a 
two-dimensional  test  case  such  as  the  driven  cavity  problem  (DCP).  The 
research  documented  in  this  report  originated  out  of  an  attempt  to  apply  his 
algorithm  to  the  DCP.  It  is  shown  herein  that  Dr.  Shoosmith’ s  method  cannot 
solve  the  DCP.  In  the  course  of  the  Investigation  additional  results,  some  of 
which  may  be  new,  were  obtained. 

The  driven  cavity  problem  is  a  standard  problem  in  computational  fluid 
dynamics.  A  National  Aeronautics  and  Space  Administration  publication  in  1975 
[Rubin]  cites  10  articles  on  the  subject  during  the  period  1965  through  1974. 
This  same  publication  documents  DCP  solutions  obtained  by  eight  research  teauns 
working  in  coordination  with  NASA.  The  computations  were  typically  performed 
on  a  CDC  6400  mainframe  computer  for  the  purpose  of  comparing  different 
approaches.  Undoubtedly  many  other  papers  on  the  DCP  have  made  their  way  into 
the  literature  since  the  writing  of  that  publication.  The  primary  advantage 
of  using  this  problem  is  that  it  is  simple  to  describe  but  the  solution  is 
fairly  complicated  and  captures  many  salient  features  of  fluid  flow  phenomena. 
Another  advantage  is  that  accurate  solutions  are  readily  available  for 
comparison. 

If  successful.  Dr.  Shoosmith’ s  proposed  method  for  solving  large  systems 
of  first-order  differential  equations  would  make  use  ot  parallel  computation 
by  distributing  eigenvector  evaluations  associated  with  an  nth-order  matrix  to 
n  parallel  processors.  The  strength  of  the  method  is  its  ability  to  perform 
the  eigenvector  computation  efficiently,  but  the  method  requires  the  matrix  tc 
have  real,  distinct  eigenvalues.  His  one-dimensional  test  case  provided  a 
solution  to  Burger’s  equation  and  the  associated  matrix  did  in  fact  have  real, 
distinct  eigenvalues.  It  is  shown  in  this  report,  however,  that  the  spectrum 
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of  the  matrix  arising  in  the  driven  cavity  problem  contains  complex 
eigenvalues  and,  unfortunately,  Dr.  Shoosmlth’s  method  cannot  provide  a 
solution.  While  Investigating  this  problem,  however,  it  was  found  that  by 
computing  the  spectrum  of  the  matrix  associated  with  a  very  coarse  grid  one 
can  gain  considerable  information  about  the  solution.  It  is  this  observation 
which  may  be  new  and  of  interest  to  practicing  computational  fluid 
dynamicists. 


Figure  1:  Geometry  of  the  Driven  Cavity 
The  Driven  Cavity  Problem 

Consider  a  rectangular  cavity  (Figure  1)  with  fixed  sides  and  bottom 
which  contains  a  viscous,  incompressible  fluid.  Along  the  top  edge  of  the 
cavity  a  flat  surface,  the  "lid,"  capable  of  translating  to  the  left  is  in 
contact  with  the  fluid.  At  time  t^  the  lid  instauitaneously  begins  movement  to 
the  left  at  a  constant  speed,  U.  The  Navier-Stokes  equations  in  the  stream 
function-vorticity  formulation  are  [Ames]: 


where  <  (vorticity)  and  \lt  (stream  function)  are  related  to  the  x  auid  y 


(1) 

(2) 
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velocity  components  (u,v)  at  every  point  within  the  cavity  by 


\i  = 


V  =  - 


C  =  V  -  u 
^  X  y 


y  X 

For  purposes  of  this  report,  the  driven  cavity  problem  is  defined  as  follows: 
If  at  time  t  =  0  the  velocity  everywhere  within  the  cavity  is  zero,  determine 
<  and  0  at  every  point  within  the  cavity  for  all  time  greater  than  t  =  0. 
Additionally,  in  this  report  the  length  and  width  of  the  cavity  are  one 
distance  unit  and  the  lid  speed  is  one  distance  unit  per  unit  time.  With 
these  values  the  Reynolds  number,  R  =  VL/v,  becomes  simply  1/v  where  v  is  the 
fluid’s  kinematic  viscosity. 


Overview  of  Report 

The  following  section  describes  how  the  partial  differential  equation  for 
vorticity,  C  .  Is  converted  into  a  system  of  ordinary  differential  equations. 
Section  3  provides  solutions  to  the  driven  cavity  problem  obtained  by  use  of 
conventional  methods.  Dr.  Shoosmith’s  method  for  solving  the  coupled  system 
of  ordinary  differential  equations  is  outlined  in  Section  4.  In  Section  5  the 
spectral  properties  of  the  coefficient  matrix  are  presented  and  discussed. 

The  initial  eigenvalues  and  their  migration  during  the  time  evolution  of  the 
solution  are  remarkable  because  they  strongly  indicate  that  complex 
eigenvalues  arise  in  ail  numerical  solutions  of  the  DCP  and  thus  the  method  of 
Dr.  Shoosmith  cauinot  be  applied.  Further  interesting  consequences  of  the 
eigenvalues  of  A  and  their  migration  as  the  solution  proceeds  to  steady  state 
are  discussed  as  well.  Section  6  contains  some  concluding  remarks  including 
suggestions  on  further  research. 


Computational  Equipment  and  Software 

All  computations  were  performed  on  a  personal  computer,  a  Unisys  Series 
800/20C  with  an  80386  CPU  and  an  80387  coprocessor.  The  programs  were  written 
using  Borland’s  Turbo  Pascal  Version  6,  and  are  available  from  the  author  upon 
request. 
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Section  2 


Applying  the  Method  of  Lines  to  the  Partial  Differential  Eqtiation 

The  method  of  lines  is  a  semi-discrete  method  which  is  sometimes  used  to 
solve  partial  differential  equations  [Ames].  The  basic  idea  is  to  discretize 
all  independent  variables  but  one,  thus  converting  the  system  of  partial 
differential  equations  into  a  system  of  ordinary  differential  equations.  The 
driven  cavity  problem  is  amenable  to  this  approach. 


X 

Figure  2:  Grid  for  Driven  Cavity  Problem 


The  method  of  lines  requires  a  discretization  of  two  of  the  three 
independent  variables  and  the  usual  procedure  is  to  select  these  to  be  the 
spatial  variables.  Consider  a  regular  grid  superimposed  on  the  interior  of 
tl.5  cavity  as  depicted  in  Figure  2.  We  denote  the  number  of  interior  points 
on  any  given  row  by  N^,  and  the  number  of  interior  points  on  any  given  column 
by  N^.  The  points  are  labeled  such  that  the  point  is  located  at 
X  =  iAx,  y  =  JAy.  Discretization  gives  rise  to  new  dependent  variables  at 
each  grid  point: 

=  C(iAx,  JAy)  =  0(iAx,  JAy) 

where 

£  =  (i  -  1)N^  +  J 
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Note  that  this  numbering  system  produces  n  =N  N  new  variables  to  represent 

X  y 

both  C  and  0.  It  is  convenient  to  slightly  abuse  the  notation  and  define 
the  column  vectors  of  these  new  unknowns  as  follows: 


c  = 

• 

• 

P  = 

r  ^] 

• 

• 

.<n. 

L'^ni 

It  is  straightforward,  now,  to  approximate  both  first  and  second-order 
derivatives  appearing  in  the  first  of  the  stream  fionction-vorticity  equations 
(Equation  1)  and  thereby  create  a  system  of  first-order  ordinary  differential 
equations  in  ^  : 

k  =  A(0)C  +  r(0)  (3) 

where  A(0)  is  an  n^^-order  matrix  and  is  a  column  vector  with  n  rows. 

The  derivation  of  each  of  these  will  be  be  discussed  below. 

Again  using  finite  differences  to  approximate  second-order  derivatives,  we 
have  in  lieu  of  Poisson's  equation  (Equation  2)  a  linear  algebraic  system  of 
equations: 

Dip  =  -<  (4) 

where  D  is  a  constant  n*'**-order  matrix  whose  specific  values  will  be  derived 
below.  It  is  clear,  now,  that  we  must  solve  a  coupled  system  of  equations, 
one  a  system  of  first-order  ordinary  differential  equations  emd  the  other  a 
linear  algebraic  system  of  equations. 


,5 
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Figure  3:  Node  Types  for  Driven  Cavity  Grid 

Determination  of  the  components  of  the  A(0)  matrix  is  straightforward  but 
a  bit  tedious  due  to  the  many  different  cases  to  consider.  There  are,  in 
fact,  nine  separate  situations  which  can  occur  as  Figure  3  depicts.  A  sample 
derivation  for  one  of  these  cases,  the  top  center,  is  now  presented  and  the 
results  of  similar  derivations  for  the  other  eight  possibilities  are  presented 
in  Appendix  A. 


Figure  4:  Top  Interior  Grid  Point 


6 


Consider  the  grid  points  around  a  general  node  on  the  top  row,  as 
depicted  in  Figure  4.  Using  central  differences  for  the  spatial  partial 
derivatives  yields  the  following  O(h^)  semi-discrete  equation  at  node  1: 


[•"s  ■  *«-N  1 

X 

^^+l  ^£-1 

4. 

'I’m  -  *t-i 

"  ^£-N 

X 

2Ay 

2Ax 

2Ax 

2Ay 

+  Jl_ 
R 


^£-1  “  *  ^l+l 


H-n 


Ax 


Ay 


We  need  to  invoke  the  boundary  conditions,  now,  to  determine  amd  Cg.  We 

let  0  =  0  on  the  entire  boundary,  but  <g  must  be  determined.  We  assume  a  no 
slip  condition  at  the  fluid-lid  interface,  so 


d0 
~  dy 


=  -1 


A  fictitious  node  at  A  is  now  introduced  so  that 


30 

dy 


a  -  =  -1 

2Ay 


Thus, 


*A  '  -  2Ay 


An  O(h^)  discrete  approximation  to  Poisson’s  equation  at  S  is 


=  -C 


Ax 


Ay 


This  yields 


<S  =  -2 


Ay^ 
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Applying  this  result  to  the  C  equation  and  performing  some  algebraic 
manipulation  we  obtain 


4. 

1 

4AxAy  RAy^ 

x 

4AxAy  RAx^ 

H-1 


-2 


1 

1 

1 

RAx^ 

RAy^ 

4AxAy  RAx^ 

2AxAy 


RAy 


The  terms  multiplying  and  would  be  elements  of  Altp) 

on  row  ^  in  columns  ^“li  and  £+1,  respectively.  The  last  term 

corresponds  to  the  element  in  row  I  of  the  r(</>)  vector. 


The  matrix  D  in  the  discretized  Poisson  equation  is  easier  to 
derive  because  all  terms  are  constant.  The  same  nine  cases  occur,  however, 
amd  each  must  be  considered  separately.  For  example,  consider  a  general  node 
in  the  extreme  right  column  as  depicted  in  Figure  5. 
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The  discretized  Poisson  equation,  then,  becomes 


Ax 


Ay 


=  -C, 


After  some  manipulation  we  have 


-h  *t-«  *  A  *1-1  * 

Ay  X  Ax 


-2 


Ax 


Ay^  J 


*l  *  -i  Vn  ' 

Ay  X 


The  equations  for  all  nine  node  types  are  summarized  in  Appendix  A. 

This  concludes  the  translation  of  the  partial  differential  equation 
into  a  system  of  ordinary  differential  equations  and  of  Poisson’ s  equation 
into  an  accompanying,  coupled  system  of  linear  algebraic  equations.  In  the 
following  section,  this  system  is  solved  by  conventional  means  on  grids  with 
varying  degrees  of  refinement. 
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Section  3 


SOLUTIONS  TO  THE  DRIVEN  CAVITY  PROBLEM  BY  CONVENTIONAL  METHODS 


The  Method 

In  this  section  results  from  solving  the  two  systems  of  equations  derived 
above  are  presented.  Those  equations,  repeated  here  for  convenience,  were 


C  =  A(0)C  +  7(0) 
D0  =  -< 


(3) 

(4) 


where  A(0)  is  a  matrix  of  order  n  (n  =  N  N  )  which  depends  upon  the  value  of 

X  y 

0,  7(0)  is  a  column  vector  with  n  components  which  also  depends  upon  the  value 
of  0,  and  D  is  a  constant  matrix  of  order  n.  Detailed  expressions  for  A(0) 
and  7(0)  can  be  found  in  Appendix  A. 


Solutions  of  Equations  3  and  4  were  obtained  by  use  of  the  classical 
fourth-order  Runge-Kutta  algorithm  and  the  Gauss-Seidel  iterative  method, 
respectively.  Initial  values  of  zero  were  used  for  ^  and  0.  A  single  step  of 
the  procedure  consisted  simply  of  propagating  forward  one  time  step  for  <  by 
Runge-Kutta,  then  solving  for  0  by  Gauss-Seidel  to  obtain  the  new  value  of  0. 
This  updated  value  was  then  used  to  modify  A(0)  and  7(0)  and  the  next  time 
step  could  then  be  taken.  A  fixed  time  step  was  used  for  each  run  of  the 
program.  A  stopping  time  was  established  by  allowing  the  process  to  run  until 
a  reasonably  small  value  of  the  max  norm  of  the  right  hand  side  of  the 
vorticity  equation  was  noted  (typically  IxlO'®).  It  should  be  noted  that  far 
more  efficient  differential  equation  solvers  amd  linear  algebraic  system 
solvers  are  available.  The  intent  was  not  to  solve  the  DCP  in  the  most 
efficient  manner,  but  merely  to  generate  solutions  which  could  be  used  to 
compute  the  spectrum  of  the  coefficient  matrix  A(0),  and  which  could  be 
compared  with  solutions  generated  by  the  method  of  Dr.  Shoosmith. 
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Results 


Validation  of  the  solution  method  was  accomplished  by  comparing  ^  ^^d  ip 
with  published  results.  Memory  limitations  of  Borland’s  Turbo  Pascal  and 
computer  time  limitations  prevented  runs  with  grids  as  fine  as  desired. 
However,  by  creating  a  series  of  three  runs  with  successively  finer  grids  it 
was  possible  to  numerically  validate  the  method.  In  Figure  6,  graphs  at  three 
y  locations  (y  =  0.25,  y  =  0.5,  and  y  =  0.75)  are  presented  for  both  ^  and  iji. 
Appendix  B  contains  tabular  listings  of  the  results  from  this  research.  The 
number  of  Intervals  for  the  grids  used  in  the  present  research  are  4,  8,  and 
16  which  correspond  to  n  =  9,  49,  and  225.  In  [Rubin]  a  grid  with  64 
intervals  is  used  to  determine  <  and  0  ,  and  values  for  x  and  y  increments  of 
1/16  of  a  unit  are  published.  These  values  are  presented  on  the  same  graphs. 
It  appears  that  as  grids  are  made  finer  the  values  of  C  and  computed  by  the 
method  described  above,  approach  the  pi^blished  values  of  Rubin.  Finer  grids 
could  of  course  be  used  in  order  to  be  more  certain  that  the  solution 
technique  in  this  research  is  accurate.  The  trend  is  sufficiently  convincing, 
however,  that  spectral  results  given  later  in  this  report  are  believed  to  be 
reasonable. 

A  conventional  technique  for  solving  the  DCP  having  been  presented,  along 
with  results  obtained  from  applying  the  technique,  we  turn  now  to  the  method 
of  Dr.  Shoosmith  for  solving  systems  of  first-order  ordinary  differential 
equations. 
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ttraam  function  ,  •beam  function  straam  function 


Section  4 


SOLUTION  BY  USE  OF  THE  MATRIX  EXPONENTIAL 

This  section  describes  the  approach  suggested  by  Dr.  Shoosmith 
[Shoosmith] .  The  method  involves  solving  Equation  3,  a  coupled  system  of 
first-order  differential  equations,  by  use  of  the  matrix  exponential.  A 
sketch  of  how  Shoosmith* s  method  might  be  used  to  solve  the  above  vector 
equations  follows: 

1.  Use  the  Initial  values  of  ^  and  C  to  construct  the  matrix  A(^)  and 
the  vector  y(0). 

2.  Determine  e*^,  the  "matrix  exponential",  which  allows  one  to  compute 
C  at  the  (k+l)-st  time  interval,  given  <  at  the  k-th  time 
interval  by  the  equation 

+  A'"(  -  I  )\ 

3.  Use  equation  (4)  to  solve  for  0 

4.  Recompute  the  matrix  A(0)  and  vector  y(0)  using  the  updated  \fi. 

5.  Iterate  steps  2.  through  4.  until  t 

final 

A  Pascal  program,  DCMEx  (Driven  Cavity  by  Matrix  Exponential) ,  was  written  to 
implement  this  solution  scheme  [Stafford].  Some  specifics  on  how  the 
algorithm  described  above  were  Implemented  are  now  presented. 

The  eigenvalues  of  A(0)  are  assumed  to  be  real  and  distinct  for  all  0. 
They  are  computed  by  converting  the  A  matrix  to  upper  Hessenberg  form  and  then 
applying  the  QR  algorithm.  The  matrix  A  is  then  factored  by  a  similarity 
transformation.  The  specific  technique  used  is  to  (1)  explicitly  compute  the 
right  and  left  eigenvectors,  (2)  scale  the  left  eigenvectors  such  that  the 
product  of  each  right  and  left  eigenvector  (for  the  same  eigenvalue)  is  unity, 
and,  (3)  form  two  matrices  S  and  T  with  the  columns  of  S  consisting  of  the 
right  eigenvectors  of  A(0)  and  the  rows  of  T  consisting  of  the  scaled  left 
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eigenvectors  of  A(^) .  This  results  in  T  being  an  approximation  to  the  inverse 
of  S,  and  the  matrix  A(^)  is  thus  diagonalized — i.e.,  A(0)  can  be  written  as 
the  product 

A  =  S  A  S"^ 

with  A  being  a  diagonal  matrix  with  the  eigenvalues  of  A  along  its  diagonal. 
The  eigenvector  computation  is  performed  using  the  inverse  power  method  (also 
called  "inverse  iteration")  and,  except  for  the  first  time  step,  the  program 
uses  the  previous  eigenvectors  and  eigenvalues  as  initial  guesses  for  the  new 
eigenvectors  and  eigenvalues.  The  number  of  iterations  needed  by  the  inverse 
iteration  procedure  is  used  to  control  the  step  size,  the  assumption  being 
that  if  a  pre-determined  number  of  iterations  is  required  then  too  large  a 
step  must  have  been  attempted. 

This,  then,  is  a  general  description  of  the  program  written  to  implement 
the  method  proposed  by  Dr.  Shoosmith.  Unforttinately,  the  program  was 
unsuccessful.  When  an  equal  number  of  row  euid  column  grid  points  were  used 
the  program  was  unable  to  factor  A(\^)  at  the  outset.  It  was  suspected,  and 
later  confirmed,  that  the  difficulty  arose  due  to  the  existence  of  repeated 
eigenvalues.  When  an  unequal  number  of  horizontal  and  vertical  grid  points 
were  used  the  program  operated  successfully  for  several  iterations  but 
eventually  bogged  down  well  before  reaching  a  steady  state.  Checks  of  the 
eigenvalues  seemed  to  indicate  that  although  they  begem  as  real  and  distinct, 
the  A(^)  values  evolved  in  time  to  bring  some  pairs  of  eigenvalues  close 
together.  As  these  pairs  came  sufficiently  close  together  they  again  caused 
problems  in  the  inverse  iteration  process  for  finding  the  eigenvectors.  In 
the  following  section  the  conventional  method  described  in  Section  3  is  used 
to  compute  the  spectrum  of  A(^)  as  the  solution  progresses  towards  steady 
state.  This  allows  a  better  understanding  of  why  Shoosmith’ s  method  fails  on 
this  problem. 
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Section  5 


SPECmUM  ANALYSIS  OF  THE  COEFFICIENT  MATRIX 

Using  the  conventional  methods  to  solve  the  driven  cavity  problem 
described  in  Section  3,  and  standard  matrix  eigenvalue  computation  techniques 
(reduction  to  upper  Hessenberg  form  and  OR  iteration),  the  spectrum  of  the 
coefficient  matrix  in  Equation  3  was  obtained  at  several  time  instances  during 
the  solution  process.  Eigenvalue  computations  were  carried  out  for  grids 
having  9,  25,  and  49  internal  points  (which  correspond  to  interval  widths  of 
1/4,  1/6,  zuid  1/8  unit,  respectively),  and  the  results  are  documented  in 
Tables  1,  2,  and  3,  and  depicted  graphically  in  Figures  7,  8,  and  9. 

Several  observations  can  be  made.  First,  the  eigenvalues  of  A(^)  at 
t  =  0  are  real  and  repeated  for  each  grid  coarseness.  Second,  the  eigenvalues 
come  off  the  real  axis  and  become  complex  immediately  after  t  =  0.  (Although 
the  first  value  of  t  displayed  in  the  tables  and  figures  is  at  t  =  1, 
computations  at  the  very  first  time  interval  in  the  simulation  show  the  same 
characteristic. )  Third,  note  that  regardless  of  the  grid  coarseness  the 
eigenvalue  with  largest  real  part  (smallest  real  part  in  absolute  value)  is 
around  -0.2  and  remains  reasonably  close  to  this  value  from  t  =  0  throxigh  t  = 
20  where  the  system  is  essentially  at  steady  state. 

The  first  two  observations  above  Justify  the  claim  in  the  previous 
section  that  the  method  of  Dr.  Shoosmlth  cannot  be  successful  in  solving  the 
driven  cavity  problem  because  the  eigenvalues  are  not  real  and  distinct  for 
all  time.  In  fact,  they  are  not  real  and  distinct  for  any  interval  during  the 
solution  if  the  grid  contains  an  equal  number  of  horizontal  and  vertical 
intervals.  Results  for  an  unequal  number  of  horizontal  aund  vertical  grid 
points  (N^=  4,  3)  are  displayed  in  Table  4.  This  table  implies  that 

although  the  spectrum  begins  with  real  and  distinct  roots,  as  time  progresses 
some  of  the  roots  coalesce  and  become  complex  conjugate  pairs.  This 
corresponds  precisely  to  the  behavior  exhibited  by  the  program  which 
implemented  Dr.  Shoosmlth’ s  method.  The  program  ran  nicely  until  the 
eigenvalues  had  nearly  coalesced  but  then  stalled  as  it  was  unable  to  compute 
eigenvectors  when  the  spectrum  contained  nearly  repeated  eigenvalues.  The 
third  observation  is  remarkable  because  of  its  implications  for  computational 
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TABLE  1:  SPECTRUM  ANALYSIS  OF  4X4  GRID 


t  =  0.0 

t  =  1.0 

t  =  20.0 

Real 

Imag 

Real 

Imag 

Real 

Imag 

-0.1875 

0.0000 

-0.1902 

0.0000 

-0.1987 

0.0000 

-0.4137 

0.0000 

-0.4176 

-0.0521 

-0.4374 

-0.1420 

-0.4137 

0.0000 

-0.4176 

0.0521 

-0.4374 

0.1420 

-0.6400 

0.0000 

-0.6400 

0.0000 

-0.6400 

0.0000 

-0.6400 

0.0000 

-0.6400 

-0.0737 

-0.6400 

-0.1785 

-0.6400 

0.0000 

-0.6400 

0.0737 

-0.6400 

0.1785 

-0.8663 

0.0000 

-0.8624 

0.0521 

-0.8426 

0.1420 

-0.8663 

0.0000 

-0.8624 

-0.0521 

-0.8426 

-0.1420 

-1 .0925 

0.0000 

-1 .0898 

0.0000 

-1.0813 

0.0000 

TABLE  2: 

SPECTRUM  ANALYSIS  OF  6X6  GRID 

t  =  0.0 

t  =  1.0 

t  =  20.0 

Real 

Imag 

Real 

Imag 

Real 

Imag 

0.1929 

0.0000 

-0.2291 

0.0000 

-0.2360 

0.0000 

0.4565 

0.0000 

-0.5049 

0.1152 

-0.5371 

0.2055 

0.4565 

0.0000 

-0.5049 

-0.1152 

-0.5371 

-0.2055 

0.7200 

0.0000 

-0.8163 

0.2185 

-0.8103 

0.0000 

0.8165 

0.0000 

-0.8163 

-0.2185 

-0.8743 

-0.3145 

0.8165 

0.0000 

-0.8548 

0.0000 

-0.8743 

0.3145 

1.0800 

0.0000 

-1.1560 

-0.3382 

-1.1468 

0.0719 

1.0800 

0.0000 

-1.1560 

0.3382 

-1.1468 

-0.0719 

1.1765 

0.0000 

-1.1874 

0.0186 

-1 .2020 

0.4338 

1.1765 

0.0000 

-1.1874 

-0.0186 

-1.2020 

-0.4338 

1.4400 

0.0000 

-1.4090 

0.0000 

-1.3289 

0.0000 

1 .4400 

0.0000 

-1.4400 

0.0000 

-1.4400 

-0.6518 

1 .4400 

0.0000 

-1 .4400 

-0.4491 

-1.4400 

0.0000 

1 .4400 

0.0000 

-1 .4400 

0.4491 

-1.4400 

0.6518 

1.4400 

0.0000 

-1.4710 

0.0000 

-1.5511 

0.0000 

1 .7035 

0.0000 

-1.6926 

0.0186 

-1.6780 

-0.4338 

1.7035 

0.0000 

-1.6926 

-0.0186 

-1.6780 

0.4338 

1 .8000 

0.0000 

-1.7240 

-0.3382 

-1.7332 

0.0719 

1.8000 

0.0000 

-1.7240 

0.3382 

-1.7332 

-0.0719 

2.0635 

0.0000 

-2.0252 

0.0000 

-2.0057 

0.3145 

2.0635 

o:oooo 

-2.0637 

-0.2185 

-2.0057 

-0.3145 

2.1600 

0.0000 

-2.0637 

0.2185 

-2.0697 

0.0000 

2.4235 

0.0000 

-2.3751 

-0.1152 

-2.3429 

-0.2055 

2.4235 

0.0000 

-2.3751 

0.1152 

-2.3429 

0.2055 

2.6871 

0.0000 

-2.6509 

0.0000 

-2.6440 

0.0000 
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TABLE  3: 

SPECTRUM 

ANALYSIS  OF  8X8  GRID 

9^ 

II 

o 

b 

1 

=  1.0 

1 

=  20.0 

Real 

Imag 

Real 

Imag 

Real 

Imag 

■0.1949 

0.0000 

-0.2594 

0.0000 

-0,2644 

0.0000 

•0.4723 

0.0000 

-0.5520 

-0.1064 

-0.6353 

-0.2320 

-0.4723 

0.0000 

-0.5520 

0.1064 

-0.6353 

0.2320 

-0.7498 

0.0000 

-0.8764 

0.2768 

-0.9208 

0.3829 

■0.8876 

0.0000 

-0.8764 

-0.2768 

-0.9208 

-0.3829 

-0.8876 

0.0000 

•0.9826 

0.0000 

-0.9646 

0.0000 

-1.1651 

0.0000 

-1.2837 

0.0000 

-1 .3349 

0.0000 

-1.1651 

0.0000 

-1.3298 

-0.4653 

-1 .3373 

0.6460 

-1.3774 

0.0000 

-1.3298 

0.4653 

-1.3373 

-0.6460 

-1.3774 

0.0000 

-1 .4494 

0.0000 

-1.4724 

0.0000 

-1.5803 

0.0000 

-1.7619 

-0.0421 

-1.7498 

-0.1456 

-1.6549 

0.0000 

-1.7619 

0.0421 

-1.7498 

0.1456 

-1.6549 

0.0000 

-1.8274 

0.5932 

-1.8366 

0.8500 

-1.8673 

0.0000 

-1.8274 

•0.5932 

-1.8366 

-0.8500 

-1.8673 

0.0000 

-1.9036 

0.0000 

-1.9276 

0.0000 

-2.0702 

0.0000 

-2.2044 

0.1222 

-2.2034 

0.2307 

-2.0702 

0.0000 

-2.2044 

-0.1222 

-2.2034 

-0.2307 

-2.1447 

0.0000 

-2.2865 

0.0000 

-2.2738 

0.0000 

-2.1447 

0.0000 

-2.2976 

0.0000 

-2.3323 

0.0786 

-2.2825 

-0.0000 

-2.3090 

0.6996 

-2.3323 

-0.0786 

-2.2825 

0.0000 

-2.3090 

-0.6996 

-2.3521 

-0.9797 

-2.5600 

0.0000 

-2.4862 

0.0000 

-2.3521 

0.9797 

-2.5600 

0.0000 

•2.5600 

0.0000 

-2.5600 

-1.4696 

•2.5600 

0.0000 

-2.5600 

-0.9745 

-2.5600 

1.4696 

•2.5600 

0.0000 

-2.5600 

0.9745 

•2.5600 

0.0000 

•2.5600 

0.0000 

•2  5600 

0.2821 

-2.5600 

0.3541 

-2.5600 

0.0000 

-2.5600 

-0.2821 

-2.5600 

-0.3541 

-2.5600 

0.0000 

-2.6318 

0.0000 

-2.7679 

-0.9797 

-2.8375 

-0.0000 

-2.8110 

-0.6996 

-2.7679 

0.9797 

-2.8375 

0.0000 

-2.8110 

0.6996 

-2.7877 

0.0786 

-2.9753 

0.0000 

-2.8224 

0.0000 

-2.7877 

-0.0786 

-2.9753 

0.0000 

-2.8335 

0.0000 

-2.8462 

0.0000 

•3.0498 

0.0000 

-2.9156 

-0.1222 

-2.9166 

0.2307 

-3.0498 

0.0000 

-2.9156 

0.1222 

-2.9166 

-0.2307 

-3.2527 

0.0000 

-3.2164 

0.0000 

-3.1924 

0.0000 

-3.2527 

0.0000 

-3.2926 

-0.5932 

-3.2834 

0.8500 

-3.4651 

0.0000 

-3.2926 

0.5932 

-3.2834 

-0.8500 

-3.4651 

0.0000 

-3.3581 

0.0421 

-3.3702 

0.1456 

-3.5397 

0.0000 

-3.3581 

-0.0421 

-3.3702 

-0.1456 

-3.7426 

0.0000 

-3.6706 

0.0000 

•3.6476 

0.0000 

-3.7426 

0.0000 

-3.7902 

0.4653 

-3.7827 

0.6460 

-3.9549 

0.0000 

-3.7902 

-0.4653 

-3.7827 

-0.6460 

•3.%49 

0.0000 

-3.8363 

0.0000 

-3.7851 

0.0000 

-4.2324 

0.0000 

-4.1374 

0.0000 

-4.1554 

0.0000 

-4.2324 

0.0000 

-4.2436 

-0.2768 

-4.1992 

•0.3829 

-4.3702 

0.0000 

-4.2436 

0.2768 

-4.1992 

0.3829 

-4.6477 

0.0000 

-4.5680 

-0.1064 

-4.4847 

-0.2320 

-4.6477 

0.0000 

-4.5680 

0.1064 

-4.4847 

0.2320 

•4.9251 

0.0000 

-4.8606 

0.0000 

-4.8556 

0.0000 
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Figure  7;  Eigenvalues  for  4x4  Grid 
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fluid  dynamicists. 


One  can  consider  the  system  of  first-order  ordinary  differential 
equations  as  a  time  varying  linear  system  with  a  time-varying  input  vector. 
This  is  reasonable  since  ip  varies  with  time,  and  since  A  and  y  vary  with  0, 
they  are  uniquely  determined  at  every  instant  of  time.  A  crude  but  perhaps 
useful  approach  to  understanding  the  implication  of  the  eigenvalue  with 
smallest  real  part  in  absolute  value  is  to  partition  the  time  interval  from 
t=0  to  the  time  when  steady  state  is  agreed  to  have  occurred,  [0,  t  ],  into 

S8 

intervals  over  which  A  is  assvimed  to  be  constant.  On  each  of  these  intervals 

(Tt 

the  solution  for  ^  will  contain  a  mode  of  the  form  v^e  where  a  has  a  value 

close  to  -0.2.  If  on  any  time  interval  the  forcing  function  y  excites  this 

mode,  it  is  reasonable  to  expect  that  steady  state  will  not  occur  until  the 
<rt 

time  function  e  has  fallen  to  some  reasonably  small  value.  All  other  modes 

contain  exponentials  with  more  negative  coefficients  than  -0.2  and  hence  their 

effects  would  die  out  sooner  than  the  dynamics  of  this  mode.  Thus,  the  effect 

of  this  particular  mode  will  tend  to  govern  the  time  to  steady  state.  Using  a 

<rt 

nominal  value  of  0.01  for  e  we  can  make  a  crude  estimate  of  the  time  for 
steady  state  to  occur  as  follows: 


(Ttss 

e 


a  0.01 


tss  “ 


2  In(lO) 
<r 


at  23  sec 


where  a  value  of  <r  =  -0.2  has  been  used.  Simulations  in  this  report  were  all 
run  to  final  times  of  20  seconds,  a  value  chosen  somewhat  arbitrarily  but 
based  loosely  on  the  norm  of  the  right  hand  side  of  Equation  3  becoming 
sufficiently  small  (approximately  lO"^).  The  value  of  23  seconds  predicted 
above  is,  then,  a  reasonable  estimate  of  t  for  this  Reynolds  number. 

SS 

Simulations  were  conducted  with  different  Reynolds  numbers  and  estimates  of 
times  to  steady  state  were  obtained  by  the  same  method  as  above.  Those 
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estimates  are  contained  in  Table  5.  Note  that  increasing  the  Reynolds  number 


Table  5:  Estimates  of  Times 

to  Steady  State 

for 

Various  Reynolds  Numbers 

Reynolds  Number 

<r 

Estimated  t 

SS 

1 

-19.49 

0.24 

10 

-  1.949 

2.4 

100 

-  0.1949 

24.0 

1000 

-  0.01949 

240.0 

by  a  factor  of  10  implies  an  increase  by  a  factor  of  10  in  the  time  to  steady 
state.  Observations  made  during  runs  confirmed  that  the  estimates  in  the 
table  were  reasonable,  though  no  detailed  numerical  experiments  were 
perf  ormed . 

It  is  not  difficult  to  show  that  the  eigenvalues  at  the  initial  time  vary 
in  an  uncomplicated  fashion  with  Reynolds  number.  Let  k(\p;K)  represent  the 
coefficient  matrix  A(0)  with  Reynolds  number  R.  By  examining  the  equations 
we  see  that  we  can  write 

A(0;R)  =  K 
RAx 

where  K  is  a  5-diagonal  matrix  with  =  -4  and  all  other  nonzero  elements 
unity.  The  spectrum  of  K  is  easy  to  determine  numerically.  For  the  3x3  grid 
we  have  the  following  eigenvalues  of  K: 

X  «  -1.172 
1 

X  =  X  -2.586 

2  3 

A  =  A  =  A  s*  -4.000 

4  5  6 

A  =  A  «  -5.414 

7  8 

A  «  -6.828 

9 

The  eigenvalues  at  t  =  0  for  any  Reynolds  nximber  can  be  determined  by 

multiplying  these  values  by  — ^ —  .  Thus,  the  smallest  eigenvalue  in  absolute 

RAx^ 

value  for  R  =  1,  10,  100,  and  1000  for  the  3x3  grid  becomes  or  =  -18.75, 

-1.875,  -0.1875  and  -0.01875.  Comparing  these  values  with  those  in  Table  5  we 
see  that  the  change  from  initial  to  steady  state  in  this  critical  value  is 
always  in  the  same  proportion. 


23 


Section  6 


CONCLUSION 

This  report  has  demonstrated  conclusively  that  Dr.  Shoosmith’s  method  for 
parallelizing  the  numerical  solution  of  systems  of  ordinary  differential 
equations  has  serious  limitations.  This  is  unfortunate,  since  there  is  still 
no  completely  satisfactory  method  for  taking  advantage  of  parallel  computers 
in  solving  systems  of  ODEs.  As  is  usual  with  research,  however,  more 
questions  have  been  raised  than  were  answered  during  this  investigation.  Some 
of  these  questions  and  suggestions  for  further  investigation  are  now  offered. 

Suggestions  for  Further  Research 

At  the  heart  of  the  Shoosmith  method  is  the  use  of  inverse  iteration  to 
solve  for  eigenvectors  of  the  coefficient  matrix,  A(^^»),  in  parallel.  It  is 
perhaps  possible  to  develop  an  efficient  algorithm  which  would  (1)  detect  the 
situation  where  eigenvalues  are  coalescing,  (2)  propagate  the  solution  by  an 
alternative  method  until  the  eigenvalues  become  complex,  (3)  compute 
eigenvectors  or  perhaps  pairs  of  eigenvectors  associated  with  these  complex 
roots  in  parallel,  and  (4)  continue  propagating  the  solution  in  parallel. 

This  appears  to  be  a  significant  challenge,  but  tackling  small  systems  (the 
3x3  grid  used  in  this  report,  for  example)  might  provide  insight  into  how  one 
could  proceed.  If  this  is  not  possible,  the  method  of  Dr.  Shoosmith  will  most 
likely  never  be  competitive  with  more  conventional  serial  methods  because  the 
classes  of  problems  it  would  solve  would  be  too  restricted. 

Other  techniques  for  parallelizing  the  solution  of  systems  of  ODEs  are 
being  actively  pursued,  most  notably  the  Kryolov  subspace  methods  proposed  in 
[Gallopoulos] ,  This  method  was  actually  applied  to  the  driven  cavity  problem 
in  [Saad],  but  in  the  primitive  variable  formulation  of  the  Navier  Stokes 
equations.  It  would  be  interesting  to  apply  the  Krylov  subspace  method  to.  the 
stream  function-vorticity  formulation  of  the  Navier-Stokes  equations.  A 
comparison  between  the  best  conventional  methods  and  a  Krylov  subspace  method 
using  the  same  grid  coarseness  might  prove  interesting. 
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The  observations  given  in  Section  S  relating  to  estimating  the  time  to 
steady  state  need  to  be  studied  on  more  complicated  geometries  than  just  that 
of  the  driven  cavity.  The  goal  of  such  research  would  be  to  determine  if  the 
time  to  steady  state  is  relatively  insensitive  to  the  shape  of  the  boundary. 
Since  computation  of  the  eigenvalues  of  the  matrix  on  a  coarse  grid  is 

computationally  Inexpensive,  one  could  decide  if  solving  for  the  steady  state 
flow  field  by  a  time  marching  algorithm  would  be  feasible  for  a  given  grid 
coarseness.  Perhaps  other  methods  for  predicting  the  time  to  steady  state  for 
various  Reynolds  numbers  already  exist,  but  if  not  the  results  of  this 
research  should  be  communicated  to  the  computational  fluid  dynamics  community. 

There  are  other  trends  to  be  observed  from  the  spectrum  euialysis 
presented  in  Section  5.  For  example,  the  largest  value  of  the  imaginary  parts 
of  eigenvalues  grew  from  runs  made  with  4x4  grids  to  6x6  grids  to  8x8  grids: 
the  values  were  approximately  0.18,  0.65,  and  1.5,  respectively.  These 
maximum  imaginary  parts  correspond  to  the  highest  frequency  component  to  be 
found  in  the  solution  for  that  particular  grid  coarseness.  This  information 
is  perhaps  useful,  but  three  points  are  not  sufficient  to  define  exactly  how 
the  highest  frequencies  grow  with  increasingly  finer  grids.  Another  trend 
could  be  explored.  It  was  noted  that  the  time  step  required  for  stability  in 
the  conventional  method  described  in  Section  3  had  to  be  considerably  smaller 
as  the  Reynolds  number  increased  on  grids  with  the  same  coarseness. 

Typically,  in  computational  fluid  dymanics,  capturing  the  flow  details  at  high 
Reynolds  numbers  takes  a  finer  grid  and  thus  more  computational  time  [Rubin]. 
Consequently,  it  is  natural  to  ask  if  there  exists  any  information  in  the 
spectrum  of  A(^)  which  would  suggest  how  to  overcome  this.  Also,  can  we 
predict  a  reasonable  time  step  size  by  performing  a  preliminary  analysis  of 
the  spectrum? 

Finally,  the  work  contained  in  this  report  could  be  useful  for  several 
purposes  not  originally  envisioned.  VThile  the  driven  cavity  problem  is 
discussed  considerably  in  the  literature,  this  report  provides  a  more  detailed 
development  of  the  equations  than  one  would  generally  find  in  refereed 
articles,  thereby  making  this  interesting  problem  more  accessible  to  faculty 
and  students.  At  USAFA  the  Department  of  Mathematica  Sciences  offers  a  course 
in  numerical  analysis.  Additionally,  the  Department  of  Aeronautics  uses 
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computational  fluid  dynamics  in  several  of  their  aerodynamics  courses.  Much 
of  the  tedious  work  needed  to  develop  and  verify  the  equations  for  the  driven 
cavity  problem  has  been  done  and  documented  herein.  Consequently,  instructors 
may  find  the  driven  cavity  problem  to  be  sufficiently  challenging  and 
interesting  to  be  used  in  conjunction  with  their  classes.  Different  Poisson 
solvers  could  be  compared,  for  example,  with  the  Gauss-Seidel  method  used  in 
this  report.  Also,  performance  of  implicit  ODE  solvers  might  be  compared  with 
that  of  the  Runge-Kutta  method.  A  better  physical  understanding  of  the 
Navier-Stokes  equations  might  also  be  gained  using  the  methods  in  this  report 
to  solve  this  problem  numerically.  From  the  numerical  results,  students  can 
then  determine,  for  example,  steady  state  velocity  and  pressure  distributions 
within  the  cavity. 
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APPENDIX  A 


Equations  Arising  from  the  Method  of  Lines 
Applied  to  the  Driven  Cavity  Problem 


In  Section  2  the  driven  cavity  probiem  was  defined  and  a  sample 
derivation  of  equations  arising  from  the  method  of  lines  was  given.  In  this 
appendix  the  complete  set  of  equations  used  in  all  simulations  is  presented. 
The  symbol  R  is  used  throughout  to  denote  the  Reynolds  number.  The  symbol  h 
is  used  to  denote  the  ratio  of  the  uniform  interval  widths  in  the  y  direction. 
Ay,  to  the  uniform  widths  in  the  x  direction.  Ax:  h  =  Ay/ Ax. 


Node  at  Lower  Left  Corner: 


Av^n  V 


4h  Ax  R 


X  •' 


X/ 
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Node  on  Bottom  Interior: 
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4h  Ax  R 


'^-l^)^^+N 

'  X 


.^(  -  2(1  ♦  ]  =  - 
h  Ax  '  •  X' 


Node  at  Lower  Right  Corner: 
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Node  on  Right  Interior: 
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h^Ax 


h-^-N " 


X' 


Node  at  Upper  Right  Corner: 
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Node  on  Top  Interior: 
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Node  on  Left  Interior: 
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APPENDIX  B 


APPROXIMATE  SOLUTIONS  TO  DRIVEN  CAVITY  PROBLEM 


zeta  at  y  = 

0.25 

X 

Rubin 

4x4  Grid 

8x8  Grid 

16x16  Grid 

0.00000 

0.06250 

-0.3786 

-0.30814 

0.12500 

-0.4857 

-0.28139 

-0.38188 

0.18750 

-0.4998 

-0.38597 

0.25000 

-0.441 1 

-0.03565 

-0.20528 

-0.34008 

0.31250 

-0.3449 

-0.27175 

0.37500 

-0.2481 

-0.12678 

-0.20586 

0.43750 

-0.1763 

-0.15785 

0.50000 

-0.1393 

-0.04945 

-0.09661 

-0.13248 

0.56250 

-0.1338 

-0.12657 

0.62500 

-0.1486 

-0.09941 

-0.13332 

0.68750 

-0.1729 

-0.14698 

0.75000 

-0.2026 

-0.07218 

-0.11650 

-0.16628 

0.81250 

-0.2429 

-0.19480 

0.87500 

-0.3038 

-0.15311 

-0.23705 

0.93750 

-0.3869 

-0.29086 

1.00000 

zeta  at  y  = 

0.5 

X 

Rubin 

4x4  Grid 

8x8  Grid 

16x16  Grid 

0.00000 

0.06250 

-2.003 

-1 .67875 

0.12500 

-1.453 

-0.68215 

-1.22305 

0.18750 

-0.4999 

-0.49712 

0.25000 

0.4174 

0.53867 

0.13994 

0.17619 

0.31250 

1.034 

0.62608 

0.37500 

1.291 

0.42900 

0.83624 

0.43750 

1.304 

0.86333 

0.50000 

1.163 

0.23901 

0.38380 

0.77608 

0.56250 

0.9431 

0.62927 

0.62500 

0.6923 

0.21847 

0.45747 

0.68750 

0.4336 

0.27544 

0.75000 

0.1633 

0.01852 

0.00636 

0.07885 

0.81250 

-0.1518 

-0.15454 

0.87500 

-0.5698 

-0.29700 

-0.45846 

0.93750 

-1.146 

-0.86118 

1 .00000 
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zeta  at  y  =  0.75 


X 

Rubin 

4x4  Grid  8x8  Grid 

16x1 6  Grid 

0.00000 

0.06250 

•3.868 

-2.99578 

0.12500 

1.649 

2.60673 

1 .75368 

0.18750 

3.887 

3.58481 

0.25000 

4.109 

2.80310  2.78029 

3.80277 

0.31250 

3.748 

3.44705 

0.37500 

3.271 

2.07229 

2.94137 

0.43750 

2.785 

2.42719 

0.50000 

2.313 

0.37904  1.31124 

1 .94826 

0.56250 

1.861 

1.51369 

0.62500 

1.424 

0.74053 

1.11720 

0.68750 

0.9865 

0.74171 

0.75000 

0.5156 

0.34140  0.26235 

0.35765 

0.81250 

-0.05295 

-0.08234 

0.87500 

-0.8273 

-0.31194 

-0.64903 

0.93750 

-1.934 

-1 .42992 

1.00000 

psi  at  y  =  0.25 

X 

Rubin 

4x4  Grid  8x8  Grid 

16x16  Grid 

0.00000 

0.06250 

4.72200E-04 

3.12544E-04 

0.12500 

2.21900E-03 

1 .68049E-03 

1.70100E-03 

0.18750 

5.24800E-03 

4.10614E-03 

0.25000 

9.07400E-03 

9.55538E-03  5.78823E-03 

7.09378E-03 

0.31250 

1.29600E-02 

1 .00698E-02 

0.37500 

1.61700E-02 

9.05850E-03 

1.24858E-02 

0.43750 

1.81900E-02 

1 .39684E-02 

0.50000 

1.87900E-02 

8.68265E-03  9.73202E-03 

1. 43601 E-02 

0.56250 

1.79700E-02 

1 .36939E-02 

0.62500 

1.59500E-02 

7.93981  E-03 

1.21371  E-02 

0.68750 

1 .30700E-02 

9.93331  E-03 

0.75000 

9.69900E-03 

3.87427E-03  4.75398E-03 

7.36243E-03 

0.81250 

6.23300E-03 

4.72794E-03 

0.87500 

3.11700E-03 

1 .55483E-03 

2.36728E-03 

0.93750 

8.59600E-04 

6.61075E-04 

1.00000 
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psi  at  y  =  0.5 

X 

Rubin 

4x4  Grid 

8x8  Grid 

0.00000 

0.06250 

3.88300E-03 

0.12500 

1.48400E-02 

1. 24731 E-02 

0.18750 

2.99400E-02 

0.25000 

4.52000E-02 

3.17669E-02 

2.85576E-02 

0.31250 

5.74600E-02 

0.37500 

6.51300E-02 

3.641 73E-02 

0.43750 

6.791  OOE-02 

0.50000 

6.62600E-02 

2.4391 8E-02 

3.49585E-02 

0.56250 

6.09900E-02 

0.62500 

5.30000E-02 

2.721 42E-02 

0.68750 

4.31400E-02 

0.75000 

3.22400E-02 

1.13257E-02 

1.65020E-02 

0.81250 

2.11600E*02 

0.87500 

1 .09700E-02 

5.95756E-03 

0.93750 

3.20300E-03 

1 .00000 

psi  at  y  =  0.75 


0.00000 

0.06250 

1.51400E-02 

0.12500 

0.18750 

4.521  OOE-02 
7.21500E-02 

4.78548E-02 

0.25000 

0.31250 

9.01300E-02 

9.98400E-02 

5.94535E-02 

6.70206E-02 

0.37500 

1.03000E-01 

6.7491 4E-02 

0.43750 

1.01300E-01 

0.50000 

0.56250 

9.56700E-02 

8.70800E-02 

3.08536E-02 

5.75064E-02 

0.62500 

0.68750 

7.61400E-02 

6.32900E-02 

4.30256E-02 

0.75000 

0.81250 

4.891  OOE-02 
3.35300E-02 

1.58792E-02 

2.67681  E-02 

0.87500 

0.93750 

1.00000 

1.82500E-02 

5.55400E-03 

1.07792E-02 
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